 % 初始化
clear
clc
warning off
tic;

load('./classification_results_frequency_5.mat')
% 并将每个频段的原始数据存储在对应cell中
num_channels = 18;
freq_bands = {1:3, 4:7, 8:12, 13:29, 30:70};
group_freq = cell(1, length(freq_bands)); % 使用cell数组存储每个频段的数据
data=combined_features(1,1:1260);
for band = 1:length(freq_bands)
    band_freqs = freq_bands{band}; % 获取当前频段的频率范围
    freq_data = zeros(num_channels, length(band_freqs)); % 为当前频段初始化存储矩阵
    
    for ch = 1:num_channels
        for f_idx = 1:length(band_freqs)
            f = band_freqs(f_idx); % 获取频率索引
            idx = (f-1) * num_channels + ch; % 索引从时域特征后的部分开始
            freq_data(ch, f_idx) = data(idx); % 将数据存入矩阵
        end
    end
    
    % 将该频段的频率数据存入对应的cell
    group_freq{band} = freq_data;
end
% 对每个频段的频率数据进行互信息计算
for band = 1:length(group_freq)
    current_row_data = group_freq{band}; % 获取当前频段的数据
    
    % 计算互信息
    mi_data = current_row_data'; % 转置数据以进行互信息计算
    num_variables = size(mi_data, 2); % 获取变量数量（列数）
    wind_size = size(mi_data, 1); % 获取样本数量（行数）
    mutual_info_matrix = zeros(num_variables, 'single'); % 初始化互信息矩阵
    
    for ii = 1:num_variables
        u1 = mi_data(:, ii); % 当前变量的列数据
        
        % 计算每一列变量与当前变量的互信息
        mutual_info_batch = arrayfun(@(col) EEG_calc_mi(u1, mi_data(:, col), wind_size), ii:size(mi_data, 2));
        
        % 将互信息值填充到互信息矩阵中
        mutual_info_matrix(ii, ii:end) = mutual_info_batch;
    end
    
    % 因为互信息矩阵是对称的，所以复制上三角矩阵到下三角
    mutual_info_matrix = mutual_info_matrix + mutual_info_matrix' - eye(num_variables) .* mutual_info_matrix;
    
    % 将下三角数据展开成向量
    %mi_vector = mutual_info_matrix(tril(true(size(mutual_info_matrix)), -1))';
    
    % 将互信息计算结果存储到对应的cell
    mi_results{band} = mutual_info_matrix;
end
save('mi_results','mi_results');
mi_lead_map(mi_results{5})
